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ABSTRACT 

We derive a description of the spatially inhomogeneous Bose-Einstein condensate which 
treats the system locally as a homogeneous system. This approach, similar to the Thomas- 
Fermi model for the inhomogeneous many-particle fermion system, is well-suited to de- 
scribe the atomic Bose-Einstein condensates that have recently been obtained experimen- 
tally through atomic trapping and cooling. In this paper, we confine our attention to the 
zero temperature case, although the treatment can be generalized to finite temperatures, as 
we shall discuss elsewhere. 

Several features of this approach, which we shall call the Thomas-Fermi-Bogolubov de- 
scription, are very attractive: 1. It is simpler than the Hartree-Fock-Bogolubov technique. 
We can obtain analytical results in the case of weakly interacting bosons for quantities such 
as the chemical potential, the local depletion, pairing, pressure and density of states. 2. The 
method provides an estimate for the error due to the inhomogeneity of the bose-condensed 
system. This error is a local quantity so that the validity of the description for a given trap 
and a given number of trapped atoms, can be tested as a function of position. We see for 
example that at the edge of the condensate, the Thomas-Fermi-Bogolubov theory always 
breaks down. 3. The Thomas-Fermi-Bogolubov description can be generalized to treat the 
statistical mechanics of the bose gas at finite temperatures. 
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1 Introduction 



The recently reported Bose-Einstein condensates of trapped neutral atoms H - represent 
the first unambiguous observations of a weakly interacting bose condensedgas. Quantita- 
tively, we can characterize the strength of the interaction by the expansion parameter in 

the perturbation treatment of the homogeneous bose gas , Vna 3 , where n is the density 
and a the scattering length of the inter-atomic potential. In the atomic-trap experiments, 

typically n ~ 10 12 — 10 14 cm~ 3 , and a ~ 1 — 5 nm, so that Vna? ~ 1CT 2 — 3 x 1CT 5 . 
Thus, in the sense of perturbation theory, the observed condensates are indeed textbook 
examples of weakly interacting systems. For the uniform bose gas, perturbation theory leads 
to simple analytical results. Although the trapped condensates can be described by means 
of the Hartree-Fock Bogolubov equations the latter approach does not lend itself to 

an analytical perturbation treatment. 

Intuitively, one expects that a many-body system whose density varies slowly in space 
can be described locally as a homogeneous system. Based on this picture, the Thomas- Fermi 
method || 7 was proposed for the calculation of the electron density in a heavy atom. Lieb 
and Simon 8 showed that the treatment is exact in the limit when the atomic number goes 
to infinity. Application to a confined Bose condensate was pioneered by Goldman, Silvera 
and Legget H, and recently reconsidered by Chou, Yang, and Yu |^0|-pl|. As pointed out 



by Kagan, Shlyapnikov, and Walraven||12||, the local-density description is valid when 

fi/hu > l , (l) 

where \x is the mean-field energy per particle, or chemical potential, and Tiuo is the zeropoint 
energy in the trap. 

In this paper, we derive such a description from first principles within the framework of 
the variational technique. We emphasize that, unlike the practice of neglecting the kinetic 
energy term in the Gross-Pitaevski equation which in the recent literature is sometimes called 
the Thomas-Fermi approximation, the resulting variational description is not limited to the 
condensate, but describes the depletion, pressure and all other thermodynamic quantities. 
Furthermore, like the uniform gas, the Thomas-Fermi theory leads to a perturbation treat- 
ment of the weakly interacting condensates, giving simple analytical expressions for these 
quantities. Another important advantage of the Thomas Fermi treatment is that it can be 
generalized to describe finite temperature systems, as we shall discuss in future work. In 
this paper we focus on the bose gas at zero temperature. 

The paper is organized as follows. In section 2, we generalize the usual Bogolubov 
transformation to describe spatially inhomogeneous condensates. In section 3, we introduce 
the Wigner representation and gradient expansion, which provide the tools to describe the 
nearly homogeneous systems and make the Thomas- Fermi approximation. The advantage of 
this systematic approach to the Thomas-Fermi approximation is that it provides an estimate 
of the error incurred by the inhomogeneity of the condensate, allowing one to estimate the 
accuracy of the Thomas Fermi results. We consider this point to be very important in view 
of the fact that some traps, depending on the potential and the number of trapped atoms, 
are too far from homogeneity to be described by a Thomas-Fermi description. In addition, 
even if the Thomas-Fermi description is valid in the middle of the trap, it breaks down at 
the edge of the condensate. In sections 4 and 5, we obtain the mean-field description of 
the bose system in the Thomas-Fermi approximation. The equations, derived within the 
framework of the variational principle, provide a fully self-consistent description, indicating 
that the Thomas-Fermi decription is by no means limited to weakly interacting systems. This 
remark can be expected to be of future importance in the light of recent experimental efforts 
to obtain condensates of higher density. Nonetheless, because of the special interest in the 
weakly interacting systems, we proceed in section 6 to derive a perturbation treatment and 
obtain analytical results for quantities such as the chemical potential, the local depletion, 
pairing and pressure. With the experimental atomic-traps in mind, we apply the results 
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of the general theory to the special case of a trapping potential that is of the type of a 
simple spherically symmetric harmonic oscillator in section 7. Finally, in section 8, we 
derive a density of states of the trapped weakly interacting condensate within the spirit of 
the Thomas-Fermi approximation. 



2 Generalized Bogolubov Transformation 

The Bogolubov quasi-particle concept [[nj provides a very elegant description of the inter- 
acting Bose-Einstein condensate. The quasi-particles are represented by creation (t)') and 
annihilation (if) operators that are linear combinations of regular single-particle creation (a)) 
and annihilation (a) operators. In treating a homogeneous system, for which we can work 
in a basis of single-particle plane-wave states of momentum k, the Bogolubov transforma- 
tion which relates the quasi-particle and regular particular operators, takes on a particularly 
simple form, 

T)l = x k a k + y k a_ k , 

Vk = x k a k + y k al k > ( 2 ) 

where for the purpose of describing the static properties of a condensate in equilibrium, we 
can limit the transformation parameters, x^, y^ to real numbers. Furthermore, the isotropy 
of the many-body system suggests that the transformation parameters only depend on the 
magnitude of the momentum, x k = £k and y\. = y k . Requiring the quasi-particle operators 

to be canonical, rj^rj^, = 5 k]k ', — vLivL' = 0; gives an additional constraint to x k 

and y k , 

x l-vl = ( 3 ) 

from which we can see that a single parameter a^, with Xk = cosher^ and y^ = sinhcrfc, 
suffices to parametrize the Bogolubov transformation (Q). In addition, with Eq.(|3]), we can 
also write the Bogolubov transformation as 

4 = x kvl - VkV--k , 

which is the inverse transformation of (0). 

It is useful to define the following quantities, the "distribution function" p and the "pair- 
ing function" A: 



(a k a k ) = - [cosh 2a k - 1] 



Pk 2 
A k = - (a k a_ k ) = - sinh2cr k , (5) 

where the brackets ( ) represent the ground state expectation value. The best values for x^ 
and y k are obtained variationally by minimizing the ground state free energy. 

As stated, the above description (l)-(3) only applies to homogeneous systems, whereas 
the treatment of a general (inhomogeneous) condensate, as we shall show below, involves a 
Bogolubov transformation that is quite different in appearance, from the homogeneous case. 
However, we can expect the results of the homogeneous treatment to describe the 'local' 
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behavior of an inhomogeneous condensate, provided the spatial variations of the condensate 
are sufficiently slow. In describing many-particle Fermion systems, this intuitive picture 
forms one of the key ingredients of the well-known Thomas-Fermi description of slowly 
varying many-particle systems. 

To arrive at a general treatment, we choose to work with boson- field operators, ^(x) 
and ^(x), an approach that offers the advantage of not having to specify a basis a priori. 
Furthermore, in the presence of a condensate, it is convenient to work with the fields ip(x) 
and V^( x )) which are displaced from the original fields ^(x) and ^(x) by the expectation 
value 0(x) of ^(x), 

tf(x) = ^(x) + </>(x), 

tft(x) = ^(x) + 0*(x), (6) 

where, for the purpose of describing the static properties of a condensate in equilibrium, 
can be taken to be real, and where ^(x) and V>'(x) are the displaced field operators which 
satisfy the canonical commutation relation, V>(x), V^(x') = <5(x — x'), and furthermore, 

$(x)> = $t( x )> = 0. (7) 

We introduce the Bogolubov transformation as a general linear transformation relating 
the displaced fields to the quasi-particle fields, £(x) and £*(x), 



^(x) = |^[x(x,z)e(z)-y(x,z)Ct( z 
^t(x) = |rf 3 4x*(x,z)f(z)-F*(x,z)e(z) 



which is the generalization of Eq.(|3|). The non-local nature of the generalized Bogolubov 
transformation @ should not be surprising — the 'homogeneous' Bogolubov transformation 
(§) can be written in the same form with the special feature, due to the homogeneity of the 
system, that X(m, y) and K(x, y) only depend on x — y. 

Requiring the quasi-particle fields to be canonical, leads to 

J d 3 z [X(x, z)X(z, y) - F(x, z)Y(z, y)] = <5(x - y) , (9) 

which is the generalization of (|3|). 

It is possible to derive equations for the inhomogeneous bose systems by variationally 
determining the best transformations X and Y, minimizing the free energy. This however, 
is not the path we choose to follow here. Instead, we manipulate the generalized Bogolubov 
transformations in a manner similar to the procedure to obtain the Wigner distribution from 
the off-diagonal single-particle density function. Once this is achieved, the steps that lead 
to a Thomas-Fermi description are known from quantum transport theory. One interesting 
new aspect of this treatment is that the central object of the theory is not a distribution 
function, which in some sense can still be regarded as an observable, but a transformation. 
Although this transformation determines the value of all observables, it is clearly not an 
observable quantity by itself. 
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3 Wigner Representation and Gradient Expansion 



Wigner showed that a quantum mechanical single-particle system, costumarily characterized 
by its wave function ^(x), can alternatively be fully characterized by a different function, 



Pw 



R, p ) = J d 3 r ^*(R + r/2)*(R - r/2) exp(ip ■ r) 



(10) 



where here — as in the rest of the paper — we work in units in which h=l. This function (pLDf) , 
known as the Wigner Distribution function, can be interpreted as a phase space distribution 
function Jill anc ^ leads to a description that is remarkably close to classical mechanics. The 
analogy with a classical phase space distribution function is not complete (for example pw 
can take on negative values), but can be justified by the fact that the quantum mechanical 
expectation value of observables are equal to the 'phase space integrals' of the corresponding 
classical quantities, weighted by (27r)~ 3 pw, 



Wl*> 



(tflpltf) 



d 3 x tf*(x)/(x)tf(x) 



(27T) 



d 3 p d 3 Rf(R) p w (R,p), 



d 3 x #*(x)p^(x) 



(27r)- 3 | d 3 p J d 3 Rpp w (R,p), 



etc. More recently, the many-particle generalization of the Wigner distribution has found 
many important applications in diverse areas such as nuclear [I3J and solid state physics 

An important motivation to work in the transformed representation of Eqs. (|lOl) , (x, x') — > 

(R,p), 



and its inverse 



A W {R, p) = J d 3 r A(R + r/2, R - r/2) exp(zp • r) 
A(x, x') = (2tt)" 3 J d 3 p AwQx- + x'] /2, p) exp(-zp ■ [x 



x'D, 



(12) 



(13) 



which shall henceforth be referred to as the Wigner representation, is that it is extraordinarily 
well suited to describe nearly homogeneous systems. This convenient feature follows from 
the gradient expansion |15[ |L6|]. The gradient expansion shows that, to first order, a 
'product' operator C(x, x') = / d 3 z A(x, z)B(z, x'), in the Wigner representation simply 
gives the algebraic product of A and B, Cvp(R, p) ~ A w (H,p)B w (Tl,p). The higher-order 
corrections to this approximation can be written as a series of terms containing successively 
higher-order derivatives in the (R, p)— coordinates, 



CV(R, p) w A W (R, p)B w (R, p) + 1 

11 3=1 



1 3 

8 i=i 



d 2 A w d 2 B w + d 2 A w d 2 B w 



dR] 



dp 2 dp 2 



dR 2 



dA w dB w dA w dB w 
dRj dpj dpj dRj 



( d 2 A w d 2 B w 
' dRj dp j dRj dp j 



+ 



(14) 



The first order correction in the gradient expansion flI4|) is {Ay/, Byy} PB , the Poisson bracket 
of A w and B w . If we know that the range of A w ana B w in p-space is of the order of p ci 
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then the magnitude of the derivatives dBw/dp and d 2 Bw/dp 2 in ([T3]) can be estimated to be 
of the order of B w /p c and B w /p 2 respectively. This approximation will allow us to obtain 
a very simple estimate of the 'inhomogeneity' error. 

At this point, we return to the generalized Bogolubov transformation, X(x, y), y(x, y), 
of the previous section. Working in the Wigner representation and expanding the 'canonicity' 
relation (|9|) between between X and Y in the manner of the gradient expansion, we find up 
to first order in the spatial derivatives, a relation that is similar to the constraint equation 
(0) of the homogeneous Bogolubov transformation, 

XjUR,p)-y£(R,p)«l. (15) 

Consequently, the general Bogolubov transform can be parametrized in the same way as the 
Bogolubov transform for the homogeneous bose gas, X w (R,,p) = cosh[er(R, p)], Y w (H,p) 
= sinh[er(R, p)], where for the slowly varying condensate, the a — parameters depend on 
the momentum and position : er(R, p). The distribution and pairing functions, p(x, x') 

= (^(x)^^')) and A(x, x') = — {ip(x.)ip(x.')) take on the following form in the Wigner 
representation: 

p w (R, p) = - [cosh(2(r(R, p)) - 1] , 

A iy (R,p) = ifrinh(2<7(R,p)). (16) 

The local a parametrization of the Bogolubov transformation is crucial to the Thomas- Fermi 
description and it is upon the validity of (O) that the Thomas-Fermi theory rests. The error 



introduced to flToj ) due to the inhomogeneity of the system can be estimated by the lowest 
order non- vanishing term in t he g radient expansion (0). Notice that the first order term 
in the gradient expansion of (p.5|) vanishes since it is the the sum of Poisson brackets of 
quantities with themselves. Consequently, the error has to be estimated from the second 
order term. 



4 Energy Density 

In the variational method, the quantity to minimize is F, the ground state free energy, which 
we can put in a 'local' form, F = J d 3 R /(R), where /(R) is the energy density. We 
achieve this result in two steps. In the first step, we shift to the Wigner representation in 
the integrand for the mean field expression for the ground state energy. In the second step, 
we notice that the short-range nature of the inter-atomic interaction renders the resulting 
integrand essentially 'local', i.e. the integrand contains only single (not double) integrals 
over the position variables. 

The ground state free energy is the expectation value of H — liN, where H is the many- 
body hamiltonian of the boson system, N, the number operator and Li, the chemical potential: 

H-fiN= J d 3 x ^(x) t /i(x)^(x) + ^Jd 3 x d 3 y ^(y) f ^(x) V(|x - y|)£(x)£(y), (17) 

where V(|x — y|) represents the inter-atomic potential and h{x) is the one-body part of the 
free energy, 

h(pc) = ~ + V eKk (x)-fx, (18) 
2m 

where V^ x t(x) is the external potential. 
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The presence of a condensate displaces the field operators \&(x) by their expectation value 
</>(x). To generate the variational free energy, we shall use the mean field approximation, 

in which terms of first and third order in ip and z/>' vanish fl7|) and the fourth order term 
factorizes as follows: 

^(y)^(x)^(xMy)> « 

A*(y, x)A(y, x) + p(y, x)p(x, y) + p(x, x)p(y, y). (19) 

The variational nature of this procedure is insured by the existence of a variational ground 
state that gives this type of factorization. In fact, the variational ground state corresponds 
to the choice of the gaussian wave functional [[J. 

The displacement of the fields and the factorization of the expectation values, although 
straightforward, gives rise to a somewhat lengthy expression for the free energy. It is then 
convenient to classify the different contributions by their order in <p and their functional 
dependence on p and A: 

1. hi is the one-body contribution of zeroth order in <j) to the ground state energy, 

hi = \j d 3 x d 3 y Mx)p(y, x)5(x - y). (20) 

2. In analogy with the Hartree-Fock theory, we call Vd; r given below, the direct energy 
contribution to the energy, 

V Air = \jd 3 x d 3 y p(y, y)p(x, x)V(|x - y|). (21) 

3. Using the same analogy to the Hartree-Fock treatment, the exchange energy, V exc h, is 
equal to 

V exch = ^Jd 3 x d 3 yp(x,y)p(y,x)V(|x-y|). (22) 

4. Standard Hartree-Fock theory does not describe pairing and the pairing energy, Vp a i r , 

= \J d 3 x d 3 y A*(y,x)A(y,x)V(|x - y|) , (23) 

is consequently absent from the Hartree-Fock expressions. 

In second order in 0, we find contributions that can be obtained from the above terms 
by replacing either A(x, y) or p(x, y) by 0(x)0(y). 

5. For example, the one-body contribution, due to the kinetic and potential energy of the 
condensate is hf, where 

hf = J d 3 x 0(x)/i(x)0(x). (24) 

6. Vfo. is the direct contribution to the interaction energy, stemming from the interaction of 
the condensate with the particles that have been 'forced' out of the condensate (depletion), 

VL = \jd 3 x A0(y)0(y)p(x,x)y(|x-yi). (25) 

7. Similarly, V^ ch is the exchange contribution of second order in 0, 

VLl = \J*x A0(y)^(x)p(y,x)y(|x-y|). (26) 
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8. We represent the pairing energy of the condensate with the particles out of the condensate 

by ^pair; 

V* at = \f<Px A0(y)0(x)A(y,x)y(|x-y|). (27) 

9. Finally, we denote the contribution of fourth order in 0, representing the interaction 
energy of the condensate with itself, by : 

V*+ = \j d 3 x d 3 y 2 (y)0 2 (x)\/(|x - y|). (28) 

With this notation, the mean-field expression for the ground state energy reads 

F = (H-fiN) 

= hi + Vdir + K xc h + Vpair + 

hf + 2Vt + 2Kth " 2< ir + V** y (29) 

where the minus sign of the Vp air term stems from the definition of A = —{tpip). 

At this point, we introduce the Wigner representation into the integrands of the above 
contributions to the mean-field expressions for the ground state free energy. The resulting 
expressions resemble the corresponding terms for the homogeneous gas, with an additional 
label R over which is integrated. For the sake of notational convenience we introduce the 
following integration symbol J R or / , which represents the usual integral over all of space, 

/ <i 3 R, if R is a position variable or (2n)~ 3 J d 3 p, if p is a momentum variable: 

J EE (2 I )- 3 /d s P , 

L s l* R - (30) 



The terms of zero order in then give 
^ In, J P 



pw(R, p), 

Kxch = / / / pw^R.pMp-p'W^p) , 

Jr. Jp Jp' 

V paiT =111 A w ,(R,p)^;(p-p / )A w ,(R,p , ) , 

JR Jp Jp' 

Vdir = J n I J J ^p^(R-r/2,p)p w (R + r/2,p , )exp(2q-r)t;(q) , (31) 

where v is the Fourier transform of the interaction potential, v(q) = J d 3 rV(r) exp(— iq ■ r). 

The terms that are of second order in <fi can be obtained by replacing one p or A by 
00. In the Wigner representation, this procedure yields expressions that are similar to the 
corresponding terms of zero order in with p w (R,p) or A W (R, p) replaced by a function 
Qw(R, p), where 

Q W (R, p) = J 0(R + r/2)0(R - r/2) exp(ip • r). (32) 
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Notice that the contributions of second order in 0, are non-local in the sense that their ex- 
pressions contain integrals over more than one position variable. Nevertheless, if we consider 
the scale on which the physical quantities vary in space, or in momentum space, it becomes 
apparent that the non-local integrals can be approximated by local expressions. We illus- 
trate this point by considering the exchange (Vf xch ) and pairing (Vp air ) energies. The key to 
obtain local expressions is to notice that Qw(R>, p) varies with respect to p on the scale of 
Ro \ where R is the size of the condensate. On the other hand, v (p — p') varies on the 
scale of Z" 1 where l r is the range of the atom- atom interaction. Typically R 3> l r so that 
Qw(R,~p) varies much more rapidly with respect to p than t>(p — p'). In fact, when p is 
large enough to make v (p — p') significantly different from v(p'), Qw(R, p) ~ 0. Thus, we 
can replace v(p — p') by v(p') in the integrands : 



The same considerations regarding the relative magnitude of the relevant length scales 
show that we can similarly simplify the expression of the 4 interaction energy, V^, and 

the direct interaction energies, Vdi r and V^ ir . The local expressions are most easily obtained 
by considering the difference in length scales before introducing the Wigner representation. 
In coordinate space, we notice that p(x, x) rs p(y,y) if |x — y| < l r . Thus, we can replace 
p(x, x) by p(y, y) in an integrand if it is accompanied by \^(|x — y|): 



To conclude this section, we summarize the results by remarking that the Wigner repre- 
sentation and the length scale considerations bring the free energy in an almost-local form. 
We need to qualify that statement because of the appearance of the Laplacian, a non-local 
operator, in the -contribution to the energy. In fact, it is the non-locality of this term 
that gives rise to a generalized Gross-Pitaevski or non-linear Schrodinger equation (NLSE). 
The resulting (almost-local) ground state free energy is F — J d 3 R /(R), where 




(33) 




(34) 



/(R) = / 



+ Kxt(R) - A* P( R > P) + ^cxch(R) + ^dir(R) + Vpair (R) 



+0(R) 



r-v 2 



+ Kxt(R) - li 0(R) + 2^ xch (R) + 2t&(R) - 2v*(R) 



2m 



+^(O)0 4 (R), 



(35) 
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where the exchange, direct and pairing energy densities are the integrands of the correspond- 
ing interaction energy contributions to the free energy : 

We X ch(R) = \ I I p(R,p)v(p-p)p w (R,p) , 

A Jp Jp' 

^pair(R) = g/ /, ^W{R,PHP-P')A W {K,P') , 

Wdir(R) = -v(0) [ I p W /(R,p)p w (R,p') , 
£ Jp Jp' 

^xch(R) = ^ 2 ( R )/ p pw(R,p)v(p), 
< ir (R) = ^0 2 (R)/ p a w (R,pMp), 

= ^ 2 ( R M°)X p^(R,p). (36) 

Notice that the free energy and free energy density are functional of A(R, p), p(R, p) and 
0(R). In the next section we determine the equilibrium values of A(R, p), p(R, p) and 0(R) 
by minimizing F [p, A, 0; p]. 



5 Self-Consistent Mean Field Theory 

In this section we derive the self-consistent mean-field equations that describe the nearly- 
uniform bose condensate at zero temperature. In the variational method, one minimizes the 
mean-field ground state free energy F[p, A,0;p]. Writing the integrands of the different 
contributions to the mean-field free energy in the Wigner representation, followed by the 
length scale arguments of the previous section showed that F [p w , A w , p] is essentially a 
local quantity. Finally, in the Thomas- Fermi limit of a nearly-homogeneous system, p^y (R, p) 
and Avy(R, p) are parametrized by a single Bogolubov transformation parameter cr(R, p) 
in the manner of Eq. (|16|) . Thus, to describe a nearly- homogeneous system, we minimize the 
Thomas-Fermi groundstate free energy, which is obtained from the mean-field free energy, 
assuming that pw and Aw are parametrized by a (TO), F [a, <fi; p] = F [pw{c), Aw(cr), 0; p]. 
We obtain the condensate wave function 0o(R) and Bogolubov parameter o"o(R, p) that 
describe the condensate by varying a and independently to get a minimum in F : 



5F 



5(j)(R) 
5F 



a=<r , 



5a(R,p) 



, (NLSE) 



0. (37) 



The variation, 5F/5<p = 0, gives the non-linear Schrodinger Equation (NLSE). The a 
variation, 5F/5a = 0, gives an equation for <7o(R, p). From p = \ [cosh(2<r) — 1] and A = 
|sinh(2cr) flT6|), we find that dp/da = sinh(2a) and dA/da = cosh(2a), so that 5F/5a = 
is equivalent to 

Now, several terms of the NLSE, as well as the functional derivatives 5F/5A and SF/Sp (|38|), 
depend on a and O so that the resulting equations have to be solved self-consistentlyTo 
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make the self-consistent nature of the equations more explicit, we consider the cx-dependent 
contributions to the functional derivatives, 8V eyLC \J5p, SV^/Sp and 5V vai \ v / '5 A, which we shall 
call the generalized potentials, 



U CXC h(R, p) 
EWR) 



pair 



: SV exch /5p w (R,p) = / v(p-p')p w (R,p'), 

■ 5V& r /5p w (R, p) =v(0) ph/(R, p'), 

Jp' 

R,p) = SV pail /SA w (R,p) =( i;(p-p')A w (R,p') 

J d' 



(39) 



where we name the generalized potentials after the respective interaction energies of which 
they are the functional derivatives, C/ exc h is the exchange potential, Udh the direct potential 
and [/p a i r the pairing potential. Writing the distribution and pairing function in the inte- 
grands of the generalized potentials in terms of 2a, we find with ( B8p that the generalized 
potentials implicitly depend on the functional derivatives of F : 



f4xch(R, p) 

^dir(R) 

C/p air (R,p) 



P) 2 



SF/Sp 



v(0) 



(6F/6p) 2 - (SF/SA) 2 
5F/5p 



^J(5F/5p) 2 - {8F/8Af 



np-p j 2 



-5F/5A 



(SF/8p) 2 - (SF/SAf 



(40) 



where it is understood that the functional derivatives in the integrands are evaluated at R 
and p'. Functional differentiation shows that the functional derivatives of F in turn depend 
on the generalized potentials, 



5F 



SA W {K, p) 

SF 
5p w (R, p) 



^pair(R, P) 



pair 
2 



(rMp) 



p 

2m 



+ Kxt(R) - a* 



+U cxch (R, p) + U diT (K, p) + 4>\K) [v(p) + v 



(41) 



Thus, equations ( |40"D and ( p|) self-consistently determine the generalized potentials. Fur- 
thermore, there is a dependence on the condensate wave function <\>. The latter has to be 
obtained from the NLSE : 



~ + *UR) 



p + U(R) + v(0)(f) 2 (K) 



:r) = o, 



where the potential £/(R), is equal to : 

U(R) = U dir (R) + U cxch (R, 0) - U paiT {R, 0) . 



(42) 



(43) 



This potential term, which stems from the interaction of the condensate with the particles 
out of the condensate, is absent in the simplest (low density limit) form of the NLSE, usually 
encountered in the literature. 
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The equations, ( f40D ([41|) ( f42|) and (0), are the full set of self-consistent mean field equa- 
tions that describe the condensate in the Thomas-Fermi approximation. The self-consistent 
equations for the homogeneous gas Jl7|], are recovered by putting Vext = and by assum- 



ing that is independent of position so that the kinetic energy contribution to the NLSE 
vanishes. Regarding the connection with the intuitive Thomas-Fermi model, we note that \x 
and Kxt in the self-consistent mean field equations always appear as /i — V ext (TV), so that it 
is natural to define a local effective chemical potential: 

/i eff (R) =/i-K xt (R)- (44) 

In fact, this is the essence of the Thomas- Fermi description : the system is described locally 
as a homogeneous system with a position dependent effective chemical potential (f44|). 

The solutions to the fully self-consistent equations determine the expectation value of all 
(static) physical observables as a function of the chemical potential p. One observable we 
can obtain in this manner is N, the number of trapped particles, 

w <"> = f = / r //( r -p) + / r ^( r )- («) 

the inversion of which yields /i(N), from which we can cast the results for the thermodynamic 
quantities in terms of the parameter that is controled or measured in the experiment — the 
number of atoms N. 



6 Low Density Limit 



The self-consistent equations, ([40] 

^ < 1 



(f41"D fl42|) and (|4^) can be solved iteratively. In the low 
we approximate the result by the expressions obtained 



(^cxch = ^dlr = ^pal = 0> where the superscript 



density regime, where \fna° 

after a single iteration, starting from o"o°' ) - ui,< rxrh — < iln . pall . 
indicates the order of the iteration). With this first guess we solve the NLSE and obtain the 
functional derivatives (|40l), 5F/5p, 5F/5A, yielding the first-order a— parameter (38), cr^\ 
and the generalized potentials (^), , U^^, U^ r . With these single iteration expressions 
we compute the expectation values of the observable quantities. 

In solving the NLSE, we shall assume that 0(R) varies slowly enough that we can also ne- 
glect the kinetic energy operator. To make the dependence on the scattering length explicit, 
we replace the potential by a pseudopotential, 



^pseudo (r) 



(46) 



where A = 4irh 2 a/m and the derivative operator is necessary to remove the divergency in 
the ground state free energy [|nj. 

Furthermore, we shall assume that 0(R) varies slowly enough that we can also neglect 
the kinetic energy operator in solving the NLSE ([4~2"1) : 



A 



(1) (R)1 =/i cfr (R) 



(47) 



where /i e fr is the effective chemical potential 

SF (D 



The functional derivatives (ETF) are 



5A 

6F& 
Sp 



-A 



^- yUcff (R) + 2A (1) (R) 
2m 1 



(48) 
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Consequently, the single iteration value for the Bogolubov transformation parameter a is 
equal to 



tanh 



2ai 1 \R,p) 



A0 2 (R) 



(p 2 /2m) - // eff (R) + 2A0 2 (R) 

/^eff(R) 

(p 2 /2m)+/i cff (R) ' 



(49) 



which can be recognized as the dilute uniform gas result if we put // e g = \i. 

The expression for the Bogolubov parameter from Eq.(^9|) is what we would have 
obtained with an effective energy density neglecting the interaction energies of the particles 
out of the condensate, V^ir, Kxch and V^. In other words, the effective ground state energy 
is F eff = / d 3 R /efr(R), where 



/ eff (R) = / ^-/x cff (R) + 2A0 2 (R) 
Jp |_|_ 2m 

-/i efr (R)0 2 (R) + ^ 4 (R). 



p w (R,p)-A0 2 (R)A w (R, p) 



(50) 



We obtain the results for the observable quantities by calculating their expectation values 

from the single iteration of Eq. fl49|) . For example, the condensate wave function is 
determined from the NLSE : 

A0 2 (R)^ cfr (R)-f/ (1) (R), (51) 
where the potential U(R) is the sum of the generalized potentials at zero momentum (|43|), 

C/«(R) = f/ c ( 2h(R, 0) + U®(R) - UQ T (R, 0), (52) 

evaluated with the single-iteration value for a. The single-iteration values for the generalized 
potentials are computed to be : 



U, 



(i) 



cxch 



^i(R,o) 

Thus, the condensate density is 



;R,0) = 17«(R) 



A 
3^ 



[/i cff (R)] 3/2 m 3 / 2 



A 



7T- 



[/ieff(R)] rri 



3 / 2 ^,3/2 



i 2 (R) « -/XeffCR) - 

A An 



[/. eff (R)] 3/2 m 3 / 2 . 



(53) 



(54) 



The total density n(R), including the correction to </> 2 (R) (54) and the local depletion, is 
equal to 



n(R) 



(R) + 



(R) + 



p(R,p) 



An' 



[//eff(R)] m 



3 / 2 ™3/2 



i/i eff (R)- A [/ieff (R)]3/V 3 / 2 , 



(55) 
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resulting in an expression for the density, n(R), in terms of the effective chemical potential 
//eff(R-)- Inverting this relation up to first order in Vna 3 , we obtain 



At eff (R) « An(R) 



1 + 



32 /n(R)a 3 



7T 



(56) 



which, for the homogeneous case, reduces to the well-known perturbation result. 

Finally, in a similar manner, we obtain the local pressure -P(R) from the expression for 

the effective free energy density (|50|), -P(R) = — (R), 



PCR) 



A0 4 (R) 



128 



(57) 



We can then replace (jf in (|57|) by its single- iteration value (54). Furthermore, replacing 
/iefr(R) in the resulting expression by (|56|) results in a local equation of state. 

The above results illustrate an important advantage of the Thomas-Fermi description — 
by neglecting the kinetic energy operator in the NLSE we recover simple analytical expres- 
sions for most quantities. These expressions are the analogues of the perturbation results 
for the dilute homogeneous bose gas. It is then of course very important to determine the 
regime and the conditions under which these results can be trusted. 

One source of error in the theory stems from neglecting the Laplacian operator in the 
NLSE. This approximation, although convenient, is not part of the Thomas-Fermi descrip- 
tion. It is always possible to calculate the condensate wave function numerically from the 
NLSE and proceed from there with the iteration of the self-consistent Thomas-Fermi equa- 
tions (|37p! Nevertheless, if we omit the Laplacian term, we can estimate the error by cal- 
culating cli the ratio of the kinetic energy term, — V 2 0/2m, and the non-linear potential 
energy in the NLSE, \<f) 3 , 



e L (R) = | - VV/2mA0 3 



-V 2 0(R)/0(R) 



A£(R) 



(5* 



where k c (R) = [8iran(R)] 1 is the inverse of the local coherence length, k c = A ( 7 1 (R). 

Another source of error, which cannot be remedied but is truly inherent to the Thomas- 
Fermi approximation, stems from the inhomogeneity of the system. This error is also more 
difficult to estimate, and one benefit of our approach is that the gradient expansion offers a 
'handle' on this quantity. Indeed, we use the lowest order non-vanishing term in the gradient 
expansion to estimate the error. This term is of second order because the first order term 
vanishes. We estimate its magnitude ([HP by replacing the partial derivatives with respect 
to the momentum variables by fc" 1 , since k c is a measure of the range in p of the observable 
at zero temperature. The relative error for the general product of two arbitrary operators A 
and B, e- x [AB], is then given by 



ei [AB] 



8fc c 2 (R) 



V 2 A 



w 



A 



V 2 B 



B 



w 



VA W ■ VB W 
A W B W 



(59) 



The validity of the Thomas- Fermi description depends on Xjy — = 1 flT5|), so that we use 
the accuracy of this equality to test the validity of the local homogeneity description. The 
expression ( O) can also be written as exp [<t(R, p)]exp [— er(R, p)] = 1, so that we choose 
A w as exp [<j(R, p)] and B w as exp [— <r(R, p)J to estimate the relative error e x . In fact, it is 
more convenient to work with exp(4a) then exp(cr), so that we compute the inhomogeneity 
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error [exp(4a) exp(— 4a)] of exp(4cr) and divide by 4 (since the relative error of f n is simply 
n x the relative error of f). In this manner, we find that 



e* [exp(cr) exp(-cr)] 



-Ci [exp(4a) exp(-4a)] 



Vexp(4cr) 



exp(4cr) 



8Ag(R) 

With the single-iteration value for the low-density condensate, 

exp[4a(R,p)] = l + ^W 
we find that the inhomogeneity error, ej(R) (|60|), is equal to 



(60) 



(61) 



e,;(R) 



1 cxt (R)A C (R) 



(p 2 /2m) + 2/i cfr (R) 



(62) 



where F ext is the force of the external potential, F ext = —W ext . As expected, the error 
is largest for p = 0, and using the p = — value, we obtain a simple position dependent 
estimate for the inhomogeneity error, ej(R), 



1 



ei (R) = -|F cxt (R)A c (R)/A0 2 (R)|< 



(63) 



where we replaced fj, e g by A0 2 . By equating this error fl53p to a chosen value, e cut 1, 
reflecting the accuracy we demand from the theory, we can determine the spatial boundary 
beyond which the Thomas Fermi Theory is less accurate than e cut . 



7 Spherically Symmetric Harmonic Oscillator Trap 

We now specialize V^ xt (R) to a harmonic oscillator potential, 

Kxt(R) = \^{R/Lf , 
where L is the size of the harmonic oscilator ground state, 

L = 



(64) 



h 



(65) 



and compute the expectation value of important quantities in the low density limit of the 
previous section. 

In zeroth order in the iteration, we recover the results of Baym and Pethick ||18|| . From 
([17]) we see that 



[//-y ext (R)]/A 



Rn 



8iraL 4 
15 



1 - (R/R f 



(66) 



where Rq is the size of the condensate, Rq = J jj^L. In zeroth order, all particles are in the 
condensate, so that N = / R 2 (R), and 



" T {—) 



and consequently, 



Rq — L 



/15aN\ 1/5 



V L I ■ 
The local coherence length, A C (R) is given by 



A,(R) 



L 2 



Rl-R 2 



(67) 



(6* 



(69) 



Before we proceed to calculate the perturbation corrections to the observables, we test 
the validity of the low density Thomas-Fermi formalism by calculating the errors. The error 
due to neglecting the Laplacian in the NLSE, e^(R) flSBQ, is easily computed with 



e L (R) 



\Ro 



[3-2QB/i? ) 2 ] 

(i - (R/R yr 



(70) 



from which we see that the laplacian can be omitted in the NLSE on condition that the size of 

the condensate is much larger than the size of the ground state, R 3> L, or (IhaNj L) 1//5 I. 
The error due to the departure of the BEC from homogeneity, e«(R) (|60|), is 



(R) = = Nr- 



L 

R 



(R/Ro) 



(1 - (R/Ro 



12^3 



(71) 



Again, notice that ej is small over most of the condensate region (R < Rq) if Rq ^> L. 

In Fie. (1) we show the density, [<f>^(R)] 2 /[<f>^{R = 0)] 2 ©, and both errors, and 
e», as a function of the distance to the middle of the trap. The curves are calculated for a 
harmonic oscillator trap of L = 1/im and an inter- atomic interaction with scattering length 
a = 5 nm. The dotted lines correspond to iV = 10 3 atoms in the trap, and the full line gives 
the results for iV = 10 6 atoms. Notice that for 10 3 particles, the Laplacian error is already 
substantial (~ 10%) in the middle of the trap. In contrast, only at R = l.SL (to be compared 
to Rq = 2AL) does the inhomogeneity error become of comparable magnitude. This indicates 
that even for as few as 1000 particles in the trap, the Thomas-Fermi description could be 
reasonably accurate for these parameters, provided one keeps the kinetic energy term in 
solving the NLSE. For 10 6 atoms, and only become of the order of 10% at R = 9.0L, 
whereas Rq = 9AL, which shows that the Thomas-Fermi description and neglecting the 
Laplacian operator are valid approximations in almost all of the condensate region. 

Under this condition, it is meaningful to calculate the perturbation corrections to the 
expectation values of the observable quantities. Including the perturbation correction, the 
local density ( |55|) is equal to 



n(R) 



R l 



8vraL 4 



(R/Ro 



1 2 ^ aR o A 

i - ^^rVi - (R/Ro) 



The number of trapped particles, N, is obtained by integrating over the density n(R) 

rRo 



N — [ n{R) = 4tt / ° dRR 2 n(R) 
Jr Jo 



(72) 



(73) 



16 




Distance (in units of L) 

Figure 1: (a) Condensate density for N = 10 3 and 10 6 ; (b) Error incurred in neglecting kinetic term in NLSE; 
(c) Error incurred in Thomas-Fermi approximation. Length scale on horizontal axis is in units of L, the extend of 
the ground state wave function. Calculations are done for L = 10~ 4 cm, scattering length a = 5 X cm. 



17 



which leads to 



N 



1 L /2/i\ 5 / 2 /2/A 



15 a \Tluj) 



24 \hujj 



(74) 



The inverse relation, /i as a function of N, can be obtained by solving for \x iteratively in the 
previous equation, wich gives up to second iteration, the following result : 



%0J 




H 


\/2 




~2 






" 60 





ft 



Similarly, we obtain the condensate density from fl54|) or the local depletion, d(R) 



-[^(R)n/UW(R 



2^/2 aR 
3tt L 2 



1 - (i2/i?o) 



(75) 
n«(R) 

(76) 



In Fig. (2), we show the local depletion as a function of position for the same parameters as 
those of Fig. (1). The local pressure is shown in Fig. (3). 



Local Depletion 






4 6 
Distance (in units of L) 



Figure 2: Depletion, defined as d(R) = [n(R) — 2 (R)] /0 2 (R), for the same systems as Fig. (1). 



To conclude this section, we repeat that the condition for the validity of the Thomas- 
Fermi description is that the size of the condensate exceeds the size of the ground state of 
the trap, Ro 3> L. An equivalent condition is that the coherence length in the middle of the 
condensate is smaller than the size of the ground state X C (R = 0) C I, or that the chemical 
potential exceeds the ground state energy, ji ^> (hu/2). These statements do not depend on 
the details of the trapping potential. Of course, the shape of the condensate, the boundary 
where the Thomas-Fermi description breaks down, and the expectation values of the local 
observables do depend on the shape of the potential. In this section, we gave the results for 
a spherically symmetric harmonic oscillator potential. For the convenience of the reader we 
tabulate several of the results up to first non-vanishing order in table I. 
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Pressure 




Distance (in units of L) 

Figure 3: Pressure for the same systems as Fig. (1). 



8 Density of States 

In the Thomas-Fermi picture, the system is locally equivalent to a uniform system. Therefore, 
there are 'local' excitations which in the low-density regime are described by the following 
energy spectrum : 

e p (R) = + ^cff(R)) 2 - /i c 2 ff (R) + /i , (77) 

which is well known from the Bogolubov treatment of the uniform case. The local dispersion 
relation ([77]) describes a phonon with position dependent sound velocity. 

To obtain the excitation of the whole system we compute the density of states using the 
formula 

g(e) = 5(e-e t ), (78) 

i 

where Y^i represents the sum over all excited states. In the spirit of the Thomas-Fermi 
approximation we take 

g(e)= [ [ 5(e-e p (R)) . (79) 
Jr Jp 

After integration over the momentum variable, we obtain 

where p e (R) is the momentum of a particle at position R with energy e. When calculating 
the remaining integral over space, we need to distinguish between spatial region (I) with 
condensate and a second region (II) without condensate, shown schematically in Fig. (4). 

It is necessary to break up the integral (|80|) over the different integration regions, because 

the dispersion relations for the excitations are different. In region (I), we use the Bogoliubov 



dp 



— i 



(80) 
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II 



I 

Distance 



R II R e 



Figure 4: Schematic representation of the region with (region I), and without condensate (region II) for a BEC 
in a harmonic trap. The condensate density is proportional to /i e fj(R,), which is a 'mirror image' of the trapping 
potential. Particles in the condensate have energy /i and a particle excited up to energy e can move into region 
(II) as far as the classical turning point R e . 



spectrum ([77]), whereas in region (II), the atoms are essentially free particles moving in the 
trap: 



e p (R) 



V 

2m 



+ V ext (R) 



The density of states is then the sum of the integrals over region (I) and (II): 



9(e) 



(e-u) / <TR 
J(l) 



e-fj] +/4r(R) -Meff(R) 



e-lA + /4r( R ) 



r I d'RJe- Kxt(R) 



(82) 



For the special case of a spherically symmetric harmonic oscillator trap, we find the 
following expression for the density of states : 



4 fi 2 
5(e) = - 



7T (kuj) 3 



{el li -I) ( dry/T 
Jo 



(e/n — l) 2 + r 2 — r 



y/ie/y. - l) 2 + r 2 



+ 2 



e/fi 



dr r 2 Je//i — r 7 



13) 



In Figs. (5) and (6) we show the density of states for the system discussed in the previous 
section, L = 1/im, a = 5 nm, N = 10 3 (Fig. (5)) and N = 10 6 (Fig. (6)). The dotted 
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Density of States for 10 3 particles 

30 I 1 1 1 1 1 1— 




5 V 10 15 20 

Energy (in units of co/2) 



Figure 5: Density of states calculated in the Thomas-Fermi approach described in the paper. The system is a 
BEC of N = 10 3 particles interacting with a scattering length a = 5 X 10~ 7 cm, in a harmonic trap with ground 
state of extend L = 10 -4 cm. 



Density of States for 1 6 particles 




Energy (in units of <o/2) 

Figure 6: Density of states for the same system as in Fig. (6), but with N = 10 6 particles. 
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lines show the result for the interacting bose gas, the full line shows the density of states of 
the ideal gas in the same trap. The density of states starts from the chemical potential //, 
consistent with (|77|), which implies that the energies are measured from the bottom of the 
potential well so that a particle of zero momentum in the condensate has energy /i. If we 
were to set out the density of states as a function of excitation energy e — fi, the density 
of states curves for the interacting BEC-systems would be shifted to the left by an amount 
/i. In contrast to the homogeneous BEC, the density of states for the interacting case, as 
a function of the excitation energy, grows faster than the density of states of the ideal gas. 
The reason is purely geometrical: the phonon has a much larger volume in coordinate space 
available (at least the volume of the condensate) than the non-interacting boson that received 
the same amount of energy and can only move near the bottom of the potential well. This 
effect outweighs the fact that the momentum space volume available to the phonon is less 
than the momentum space volume available to the non-interacting particle with the same 
energy. 

Of course, the sharpness of the boundary between region (I) and (II), is an artifact of 
neglecting the Laplacian operator in the NLSE. Nevertheless, except for a region near the 
boundary, we argue that the rest of space is well-described and that the contribution of 
the near-boundary region is comparatively small so that the error that is introduced in 
the integral (|80D is small provided the Thomas-Fermi description is valid in most of the 
condensate region. 
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Table I 

Results for the spherically symmetric harmonic 
oscillator trap 



Size of the Condensate 


r = l(^ l n ) 1/5 


Chemical Potential 


r~ 2 \ L J 


Condensate Density 




Local Coherence Length 


X C (R)- J£— 


Local Depletion 


d(R) = 2 -£^l-(R/R ) 2 


Error due to neglecting the Laplacian 


e L {rt) - y R j (1 _ (R/Ro) 2 )a 


Error due the inhomogeneity 


r(P \_l(L\ 4 (R/Ro) 2 
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